Section: Scientific Foundations

Identification and approximation

Identification typically consists in approximating experimental data by the prediction of a model belonging to some model class. It consists therefore of two steps, namely the choice of a suitable model class and the determination of a model in the class that fits best with the data. The ability to solve this approximation problem, often non-trivial and ill-posed, impinges on the effectiveness of a method.

Particular attention is payed within the team to the class of stable linear time-invariant systems, in particular resonant ones, and in isotropically diffusive systems, with techniques that dwell on functional and harmonic analysis. In fact one often restricts to a smaller class—e.g., rational models of suitable degree (resonant systems, see section  4.2 ) or other structural constraints—and this leads us to split the identification problem in two consecutive steps:

1. Seek a stable but infinite (numerically: high) dimensional model to fit the data. Mathematically speaking, this step consists in reconstructing a function analytic in the right half-plane or in the unit disk (the transfer function), from its values on an interval of the imaginary axis or of the unit circle (the band-width). We will embed this classical ill-posed issue (i.e., the inverse Cauchy problem for the Laplace equation) into a family of well-posed extremal problems, that may be viewed as a regularization scheme of Tikhonov-type. These problems are infinite-dimensional but convex (see section 3.1.1 ).

2. Approximate the above model by a lower order one reflecting further known properties of the physical system. This step aims at reducing the complexity while bringing physical significance to the design parameters. It typically consists of a rational or meromorphic approximation procedure with prescribed number of poles in certain classes of analytic functions. Rational approximation in the complex domain is a classical but difficult non-convex problem, for which few effective methods exist. In relation to system theory, two specific difficulties superimpose on the classical situation, namely one must control the region where the poles of the approximants lie in order to ensure the stability of the model, and one has to handle matrix-valued functions when the system has several inputs and outputs, in which case the number of poles must be replaced by the McMillan degree (see section 3.1.2 ).

When identifying elliptic (Laplace, conjugate-Beltrami) partial differential equations from boundary data, point 1. above can be recast as an inverse boundary-value problem with (overdetermined Dirichlet-Neumann) data on part of the boundary of a plane domain (recover a function, analytic in a domain, from incomplete boundary data). As such, it arises naturally in higher dimensions when analytic functions get replaced by gradients of harmonic functions (see section  4.1 ). Initial motivations of the team include:

• free boundary problems in plasma control;

• the recovery of sources, that arises for instance in magneto/electro-encephalography;

• the detection of cracks and occlusions in non-destructive control.

We aim at generalizing this approach to the conjugate-Beltrami equation in dimension 2 (section 6.2.3 ) and to the Laplace equation in dimension 3 (section 6.2.1 ).

Step 2. above, i.e., meromorphic approximation with prescribed number of poles—is used to approach other inverse problems beyond harmonic identification. In fact, the way the singularities of the approximant (i.e., its poles) relate to the singularities of the approximated function is an all-pervasive theme in approximation theory: for appropriate classes of functions, the location of the poles of the approximant can be used as an estimator of the singularities of the approximated function (see section 6.2.2 ).

We provide further details on the two steps mentioned above in the sub-paragraphs to come.

Analytic approximation of incomplete boundary data

Participants : Laurent Baratchart, Slah Chaabi, Sylvain Chevillard, Yannick Fischer [Until November] , Juliette Leblond, Jean-Paul Marmorat, Jonathan Partington, Elodie Pozzi [Since October] , Fabien Seyfert.

Given a planar domain $D$, the problem is to recover an analytic function from its values on a subset of the boundary of $D$. It is convenient to normalize $D$ and apply in each particular case a conformal transformation to meet a “normalized” domain. In the simply connected case, which is that of the half-plane, we fix $D$ to be the unit disk, so that its boundary is the unit circle $T$. We denote by $H^p$ the Hardy space of exponent $p$ which is the closure of polynomials in the $L^p$-norm on the circle if $1\le p<\infty$ and the space of bounded holomorphic functions in $D$ if $p=\infty$. Functions in $H^p$ have well-defined boundary values in $L^p\left(T\right)$, which makes it possible to speak of (traces of) analytic functions on the boundary.

A standard extremal problem on the disk is [61] :

($P_0$)  Let $1\le p\le \infty$ and $f\in L^p\left(T\right)$; find a function $g\in H^p$ such that $g-f$ is of minimal norm in $L^p\left(T\right)$.

When seeking an analytic function in $D$ which approximately matches some measured values $f$ on a sub-arc $K$ of $T$, the following generalization of ($P_0$) naturally arises:

($P$)  Let $1\le p\le \infty$, $K$ a sub-arc of $T$, $f\in L^p\left(K\right)$, $\psi \in L^p(T\setminus K)$ and $M>0$; find a function $g\in H^p$ such that ${\parallel g-\psi \parallel }_{L^p(T\setminus K)}\le M$ and $g-f$ is of minimal norm in $L^p\left(K\right)$ under this constraint.

Here $\psi$ is a reference behavior capsulizing the expected behavior of the model off $K$, while $M$ is the admissible error with respect to this expectation. The value of $p$ reflects the type of stability sought after and how much one wants to smoothen the data.

To fix terminology we generically refer to ($P$) as a bounded extremal problem. The solution to this convex infinite-dimensional optimization problem can be obtained upon iteratively solving spectral equations for appropriate Hankel and Toeplitz operators, that involve a Lagrange parameter, and whose right hand-side is given by the solution to ($P_0$) for some weighted concatenation of $f$ and $\psi$. Constructive aspects are described in [45] , [47] , [12] , for $p=2$, $p=\infty$, and $1<p<\infty$, while the case $p=1$ is essentially open.

Various modifications of $\left(P\right)$ have been studied in order to meet specific needs. For instance when dealing with loss-less transfer functions (see section 4.2 ), one may want to express the constraint on $T\setminus K$ in a pointwise manner: $|g-\psi |\le M$ a.e. on $T\setminus K$, see [48] for $p=2$ and $\psi =0$.

The above-mentioned problems can be stated on an annular geometry rather than on a disk. For $p=2$ the solution proceeds much along the same lines [71] . When $K$ is the outer boundary, ($P$) regularizes a classical inverse problem occurring in nondestructive control, namely recovering a harmonic function on the inner boundary from overdetermined Dirichlet-Neumann data on the outer boundary (see sections  4.1 and 6.2 ). Interestingly perhaps, it becomes a tool to approach Bernoulli type problems for the Laplacian, where overdetermined observations are made on the outer boundary and we seek the inner boundary knowing it is a level curve of the flux (see section  6.2.3 ). Here, the Lagrange parameter indicates which deformation should be applied on the inner contour in order to improve data fitting.

Continuing effort is currently payed by the team to carry over bounded extremal problems and their solution to more general settings.

Such generalizations are twofold: on the one hand Apics considers 2-D diffusion equations with variable (but for now isotropic) conductivity, on the other hand it investigates the ordinary Laplacian in ${𝐑}^n$. The targeted applications are the determination of free boundaries in plasma control and the source detection in electro/magneto-encephalography (EEG/MEG), as well as discretization issues of the gravitational potential in geophysics (see section 6.2.2 ).

An isotropic diffusion equation in dimension 2 can be recast as a so-called conjugate or real Beltrami equation [70] . This way analytic functions get replaced by “generalized” ones in problems ($P_0$) and ($P$). Hardy spaces of solutions, which are more general than Sobolev ones and allow one to handle $L^p$ boundary conditions, have been introduced when $1<p<\infty$ [8] , [43] . The expansions of solutions needed to constructively handle such problems have been studied in [17] , [21] . The goal is to solve the analog of ($P$) in this context to approach Bernoulli-type problems (see section 6.2.1 , [22] ).

At present, bounded extremal problems for the $n$-D Laplacian are considered on half-spaces or balls. Following [82] , Hardy spaces are defined as gradients of harmonic functions satisfying $L^p$ growth conditions on inner hyperplanes or spheres. From the constructive viewpoint, when $p=2$, spherical harmonics offer a reasonable substitute to Fourier expansions [2] . Only very recently were we able to define operators of Hankel type whose singular values connect to the solution of ($P_0$) in BMO norms. The $L^p$ problem also makes contact with some nonlinear PDE's, namely to the $p$-Laplacian. The goal here is to solve the analog of ($P$) on spherical shells to approach inverse diffusion problems across a conductor layer.

Meromorphic and rational approximation

Participants : Laurent Baratchart, José Grimm, Martine Olivi, Edward Saff, Herbert Stahl [TFH Berlin] , Maxim Yattselev.

Let as before $D$ designate the unit disk, and $T$ the unit circle. We further put $R_N$ for the set of rational functions with at most $N$ poles in $D$, which allows us to define the meromorphic functions in $L^p\left(T\right)$ as the traces of functions in $H^p+R_N$.

A natural generalization of problem ($P_0$) is

($P_N$)  Let $1\le p\le \infty$, $N\ge 0$ an integer, and $f\in L^p\left(T\right)$; find a function $g_N\in H^p+R_N$ such that $g_N-f$ is of minimal norm in $L^p\left(T\right)$.

Problem ($P_N$) aims, on the one hand, at solving inverse potential problems from overdetermined Dirichlet-Neumann data, namely to recover approximate solutions of the inhomogeneous Laplace equation $\Delta u=\mu$, with $\mu$ some (unknown) distribution, which will be discretized by the process as a linear combination of $N$ Dirac masses. On the other hand, it is used to perform the second step of the identification scheme described in section 3.1 , namely rational approximation with a prescribed number of poles to a function analytic in the right half-plane, when one maps the latter conformally to the complement of $D$ and solve ($P_N$) for the transformed function on $T$.

Only for $p=\infty$ and continuous $f$ it is known how to solve ($P_N$) in closed form. The unique solution is given by AAK theory (named after Adamjan, Arov and Krein), that connects the spectral decomposition of Hankel operators with best approximation in Hankel norm  [78] . This theory allows one to express $g_N$ in terms of the singular vectors of the Hankel operator with symbol $f$. The continuity of $g_N$ as a function of $f$ only holds for stronger norms than uniform.

The case $p=2$ is of special importance. In particular when $f\in {\overline{H}}^2$, the Hardy space of exponent 2 of the complement of $D$ in the complex plane (by definition, $h\left(z\right)$ belongs to ${\overline{H}}^p$ if, and only if $h(1/z)$ belongs to $H^p$), then ($P_N$) reduces to rational approximation. Moreover, it turns out that the associated solution $g_N\in R_N$ has no pole outside $D$, hence it is a stable rational approximant to $f$. However, in contrast with the situation when $p=\infty$, this approximant may not be unique.

The former Miaou project (predecessor of Apics) has designed an adapted steepest-descent algorithm for the case $p=2$ whose convergence to a local minimum is guaranteed; until now it seems to be the only procedure meeting this property. Roughly speaking, it is a gradient algorithm that proceeds recursively with respect to the order $N$ of the approximant, in a compact region of the parameter space [41] . Although it has proved rather effective in all applications carried out so far (see sections 4.1 , 4.2 ), it is not known whether the absolute minimum can always be obtained by choosing initial conditions corresponding to critical points of lower degree (as done by the Endymion software, section  5.5 , and RARL2 software, section 5.2 ).

In order to establish convergence results of the algorithm to the global minimum, Apics has undergone a long-haul study of the number and nature of critical points, in which tools from differential topology and operator theory team up with classical approximation theory. The main discovery is that the nature of the critical points (e.g., local minima, saddles...) depends on the decrease of the interpolation error to $f$ as $N$ increases [49] . Based on this, sufficient conditions have been developed for a local minimum to be unique. These conditions are hard to use in practice because they require to establish strong estimates for the approximation error. Such estimates are often difficult to obtain for a given function, and are usually valid only for large $N$. Examples where uniqueness or asymptotic uniqueness has been proved this way include transfer functions of relaxation systems (i.e., Markov functions) [51] , the exponential function, and meromorphic functions [11] . The case where $f$ is the Cauchy integral on a hyperbolic geodesic arc of a Dini-continuous function which does not vanish “too much” has been recently answered in the positive. An analog to AAK theory has been carried out for $2\le p<\infty$ [12] . Although not computationally as powerful, it has better continuity properties and stresses a continuous link between rational approximation in $H^2$ and meromorphic approximation in the uniform norm, allowing one to use, in either context, techniques available from the other (When $1\le p<2$, problem ($P_N$) is still fairly open.).

A common feature to all these problems is that critical point equations express non-Hermitian orthogonality relations for the denominator of the approximant. This is used in an essential manner to assess the behavior of the poles of the approximants to functions with branched singularities, which is of particular interest for inverse source problems (cf. section 6.2.2 ).

In higher dimensions, the analog of problem ($P_N$) is the approximation of a vector field with gradients of potentials generated by $N$ point masses instead of meromorphic functions. The issue is by no means understood at present, and is a major endeavor of future research problems.

Certain constrained rational approximation problems, of special interest in identification and design of passive systems, arise when putting additional requirements on the approximant, for instance that it should be smaller than 1 in modulus. Such questions have become over years an increasingly significant part of the team's activity (see section  4.2 ). When translated over to the circle, a prototypical formulation consists in approximating the modulus of a given function by the modulus of a rational function of degree $n$. When $p=2$ this problem can be reduced to a series of standard rational approximation problems, but usually one needs to solve it for $p=\infty$. The case where $\left|f\right|$ is a piecewise constant function with values 0 and 1 can also be approached via classical Zolotarev problems [81] , that can be solved more or less explicitly when the pass-band consists of a single arc. A constructive solution, in the case where $\left|f\right|$ is a piecewise constant function with values 0 and 1 on several arcs (multiband filters), is one recent achievement of the team. Though the modulus of the response is the first concern in filter design, the variation of the phase must nevertheless remain under control to avoid unacceptable distortion of the signal. This is an important issue, currently under investigation within the team under contract with the CNES. From the point of view of design, rational approximants are indeed useful only if they can be translated into physical parameter values for the device to be built. This is where system theory enters the scene, as the correspondence between the frequency response (i.e., the transfer-function) and the linear differential equations that generate this response (i.e., the state-space representation), which is the object of the so-called realization process. Since filters have to be considered as dual-mode cavities, the realization issue must indeed be tackled in a $2\times 2$ matrix-valued context that adds to the complexity. A fair share of the team's research in this direction is concerned with finding realizations meeting certain constraints (imposed by the technology in use) for a transfer-function that was obtained with the above-described techniques (see section  6.6 ).

Behavior of poles of meromorphic approximants and inverse problems for the Laplacian

Participants : Laurent Baratchart, Herbert Stahl [TFH Berlin] , Maxim Yattselev.

We refer here to the behavior of the poles of best meromorphic approximants, in the $L^p$-sense on a closed curve, to functions $f$ defined as Cauchy integrals of complex measures whose support lies inside the curve. If one normalizes the contour to be the unit circle $T$, we are back to the framework of section 3.1.2 and to problem ($P_N$); the invariance of the problem under conformal mapping was established in [9] . The research so far has focused on functions whose singular set inside the contour is zero or one-dimensional.

Generally speaking, the behavior of poles is particularly important in meromorphic approximation to obtain error rates as the degree goes large and also to tackle constructive issues like uniqueness. However, the original motivation of Apics is to consider this issue in connection with the approximation of the solution to a Dirichlet-Neumann problem, so as to extract information on the singularities. The general theme is thus how do the singularities of the approximant reflect those of the approximated function? The approach to inverse problem for the 2-D Laplacian that we outline here is attractive when the singularities are zero- or one-dimensional (see section  4.1 ). It can be used as a computationally cheap preliminary step to obtain the initial guess of a more precise but heavier numerical optimization.

As regards crack detection or source recovery, the approach in question boils down to analyzing the behavior of best meromorphic approximants of a function with branch points. We were able to prove ([6] , [9] ) that the poles of the approximants accumulate in a neighborhood of the geodesic hyperbolic arc that links the endpoints of the crack, or the sources [46] . Moreover, the asymptotic density of the poles turns out to be the equilibrium distribution on the geodesic arc of the Green potential and this distribution puts heavy charge at the end points, that are thus well localized if one is able to compute sufficiently many zeros (this is where the method could fail). The case of more general cracks, as well as situations when three or more sources, require handling a finite but arbitrary number of branch points. These are outstanding open questions for applications to inverse problems (see section  6.2 ), as for the problem of a general singularity, that may be two dimensional.

Results of this type open new perspectives in non-destructive control, in that they link issues of current interest in approximation theory (the behavior of zeroes of non-Hermitian orthogonal polynomials) to some classical inverse problems for which a dual approach is thereby proposed: to approximate the boundary conditions by true solutions of the equations, rather than approximating (by discretization) the equation itself.

We wish to point out that rational or meromorphic approximation to the Cauchy transform of measure can be viewed as discretization of the logarithmic potential of that measure. If approximation takes place in the $L^p$ sense on the boundary of a domain, the discretization proceeds according to a homogeneous Sobolev $p$-norm. This formulation can be generalized to higher dimensions, even if the computational power of complex analysis is then no longer available. This makes for a long-term research project with a wide range of applications. It is interesting to mention that the case of sources in dimension three in a spherical or ellipsoidal geometry can be attacked with the above 2-D techniques, as applied to planar sections (see section  6.2 ).

Matrix-valued rational approximation

Participants : Laurent Baratchart, Martine Olivi, José Grimm, Jean-Paul Marmorat, Bernard Hanzon, Ralf Peeters [Univ. Maastricht] .

Matrix-valued approximation is necessary for handling systems with several inputs and outputs, and it generates substantial additional difficulties with respect to scalar approximation, theoretically as well as algorithmically. In the matrix case, the McMillan degree (i.e., the degree of a minimal realization in the System-Theoretic sense) generalizes the degree.

The problem we want to consider reads: Let $ℱ\in {\left(H^2\right)}^{m\times l}$ and $n$ an integer; find a rational matrix of size $m\times l$ without poles in the unit disk and of McMillan degree at most $n$ which is nearest possible to $ℱ$ in ${\left(H^2\right)}^{m\times l}$. Here the $L^2$ norm of a matrix is the square root of the sum of the squares of the norms of its entries.

The approximation algorithm designed in the scalar case generalizes to the matrix-valued situation [60] . The first difficulty here consists in the parametrization of transfer matrices of given McMillan degree $n$, and the inner matrices (i.e., matrix-valued functions that are analytic in the unit disk and unitary on the circle) of degree $n$. The latter enter the picture in an essential manner as they play the role of the denominator in a fractional representation of transfer matrices (using the so-called Douglas-Shapiro-Shields factorization).

The set of inner matrices of given degree has the structure of a smooth manifold that allows one to use differential tools as in the scalar case. In practice, one has to produce an atlas of charts (parametrization valid in a neighborhood of a point), and we must handle changes of charts in the course of the algorithm. Such parametrization can be obtained from interpolation theory and Schur type algorithms, the parameters being interpolation vectors or matrices ( [35] , [15] , [16] ). Some of these parametrizations have a particular interest for computation of realizations ([15] , [16] ), involved in the estimation of physical quantities for the synthesis of resonant filters. Two rational approximation software codes (see sections 5.2 and 5.5 ) have been developed in the team.

Problems relative to multiple local minima naturally arise in the matrix-valued case as well, but deriving criteria that guarantee uniqueness is even more difficult than in the scalar case. The already investigated case of rational functions of the sought degree (the consistency problem) was solved using rather heavy machinery [10] . The case of matrix-valued Markov functions, the first example beyond rational functions, has undergone progress only recently [40] .

Let us stress that the algorithms mentioned above are first to handle rational approximation in the matrix case in a way that converges to local minima, while meeting stability constraints on the approximant.